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Analysis  of  Different  Approaches  to  Modeling  of  Nozzle  Flows  in 

the  Near  Continuum  Regime 

E.  V.  Titov  * * * §  and  D.A.  Levin'*' 

Pennsylvania  State  University,  University  Park,  PA  16802 
N.  E.  Gimelsheiid  and  S.  F.  Girrielsheird 
ERC,  Inc.,  Edwards  AFB,  CA  93524 

Conical  nozzle  flows  are  studied  for  Reynolds  numbers  of  1,230  and  12,300  using  different 
numerical  techniques:  the  direct  simulation  Monte  Carlo  method,  Navier-Stokes/CFD  that 
accounts  for  wall  temperature  jump  and  velocity  slip,  and  statistical  and  deterministic 
approaches  of  the  BGK  equation.  Detailed  comparisons  of  the  efficiency,  stability,  accuracy, 
and  convergence  of  the  employed  numerical  techniques  provides  better  understanding  of 
their  benefits  and  deficiencies,  and  assists  in  selecting  the  most  appropriate  technique  for  a 
particular  nozzle  and  flow  application.  The  deterministic  solution  of  the  BGK  equation  was 
found  to  be  in  good  agreement  with  the  benchmark  DSMC  results,  while  there  were  some 
differences  observed  between  the  statistical  BGK  and  DSMC.  The  Navier-Stokes  solution 
differs  from  DSMC  in  the  boundary  layer. 

I.  Introduction 

Multiparameteric  analysis  of  low  and  moderate  Reynolds  number  nozzle  flows  is  often  problematic  due  to  a 
large  number  of  cases  that  need  to  be  examined  for  the  optimization  of  the  nozzle-based  device  performance. 
Experimental  studies  under  such  conditions  are  rare,  expensive,  and  often  do  not  provide  the  necessary 
precision  to  adequately  assess  nozzle  performance-related  properties  such  as  thrust,  flow  rate,  and  specific 
impulse.  The  problem  is  even  worse  for  microscale,  MEMS  fabricated  thrusters,  whose  interaction  with 
critical  spacecraft  surfaces  at  high  altitudes  needs  to  be  evaluated.  Significant  back  flow  generated  by  such 
devices  plays  a  major  role  in  the  contamination  of  optical  instruments,  solar  panels,  etc.  The  back  flow 
formation  is  very  sensitive  to  the  flow  conditions  at  the  nozzle  lip,  and  is  known  to  be  difficult  to  analyze 
experimentally. 

The  development  of  accurate  numerical  tools  capable  of  handling  micronozzle  flows  is  therefore  impor¬ 
tant,  although  it  is  complicated  by  the  change  in  the  flow  regime  from  continuum,  near  the  nozzle  throat, 
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to  transitional  closer  to  the  nozzle  exit  plane.  Both  kinetic,  such  as  the  direct  simulation  Monte  Carlo 
(DSMC)  method,  and  continuum,  mostly  based  on  the  solution  of  Navier-Stokes  (N-S)  equations,  tech¬ 
niques  meet  computational  and  physical  challenges  when  applied  to  these  flows.  The  major  problem  with 
the  DSMC  method1  is  its  often  prohibitive  computational  cost  when  high  density  portions  of  the  flow  have 
to  be  accurately  modeled.  Conventional  continuum  CFD  techniques  may  be  inapplicable  in  the  regions  of 
high  gradients  and  strong  rarefaction  even  when  the  velocity  slip  and  temperature  jump  corrections  to  the 
boundary  conditions  at  the  nozzle  surface  are  used. 

A  combined  Navier-Stokes/DSMC  approach  is  often  used  in  the  recent  years  to  model  nozzle  flows,2,3 
where  the  solution  of  the  Navier-Stokes  equations  is  used  to  model  high  density  portions  of  the  flow  inside  the 
nozzle,  while  the  DSMC  method  is  used  to  predict  more  rarefied  plume-atmosphere  interaction.  A  starting 
surface  is  used  to  hand  the  properties  obtained  in  the  continuum  modeling  off  to  the  DSMC  solver.  The  flow 
is  always  assumed  to  be  in  steady  state,  and  it  is  implied  that  the  more  rarefied,  DSMC  portions  of  the  flow 
have  negligible  effect  on  the  upstream  continuum  portions.  This  approach  becomes  however  not  applicable 
when  the  latter  implication  is  not  true  (for  example,  when  parts  of  the  hand-off  surface  are  subsonic,  or  the 
downstream  flow  impacts  the  upstream  flow  through  radiation). 

In  an  attempt  to  overcome  this  problem,  it  is  possible  to  use  a  two-way  coupled  hybrid  continuum-kinetic 
approach  (see,  for  example,  Refs.  4-6),  although,  in  addition  of  the  usual  hybridization  problems,7  the 
extantion  to  real  gases  effect  or  two-phase  models  may  not  be  straightforward.  The  deficiencies  of  the  hand- 
off  NS-DSMC  approach  become  even  more  severe  when  transient  nozzle  and  plume  flows  need  to  be  analyzed. 
In  this  case  such  an  approach,  as  well  as  two-way  coupled  hybrid  techniques,  are  difficult  to  implement.  The 
main  difficulty  is  related  to  the  temporal  changes  in  gas  properties,  most  importantly  gas  mean  free  path 
and  Mach  number,  that  necessitate  flexible,  transient  hand-off  surfaces  and  interface  boundaries. 

It  would  therefore  be  highly  desirable  to  have  a  single  method  that  allows  an  accurate  and  efficient 
one-step  modeling  of  high  density  nozzle  and  low  density  plume  flows.  One  of  the  first  attempts  to  tackle 
this  problem  was  undertaken  in  Ref.  8  where  a  statistical,  particle  approach  to  the  solution  of  the  ES-BGK 
equation  was  obtained  for  a  flow  expanding  through  a  conical  nozzle.  Generally,  the  use  of  a  simplified 
form  of  the  Boltzmann  equation,  usually  called  a  model  kinetic  equation,  such  as  Bhatnagar-Gross-Krook 
(BGK)  or  ellipsoidal  statistical  (ES-BGK)  equations,  may  be  beneficial.  A  method  based  on  the  solution  of 
model  kinetic  equations  is  expected  to  be  more  efficient  than  the  DSMC  method  in  the  continuum  and  near¬ 
continuum  regime,  and  more  accurate  than  the  solution  of  the  Navier-Stokes  equations  in  the  transitional 
regime. 

In  recent  years  there  has  been  a  renewed  interest  and  significant  advances  in  the  solution  of  model 
kinetic  equations,  with  deterministic,  either  finite  difference  or  finite  volume,  approaches  typically  used  in 
the  solution  procedure.  A  particle  approach  to  obtain  the  solution  to  the  BGK  equation  was  first  proposed 
in  Ref.  9.  It  was  then  extended  to  model  the  ES-BGK  equation  in  Ref.  10,  and  further  extended  to  include 
rotational  degrees  of  freedom  in  Ref.  8.  The  main  goal  of  this  paper  is  to  use  statistical  and  deterministic 
approaches  to  model  kinetic  equations  to  simulate  rarefied  gas  flow  through  a  conical  nozzle,  and  compare 
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the  results  in  terms  of  accuracy  and  efficiency  to  the  DSMC  method  and  the  solution  of  the  Navier-Stokes 
equations. 


II.  Geometry  and  flow  conditions 

Gas  flow  of  pure  argon  through  a  conical  nozzle  into  vacuum  is  considered  in  this  work.  A  diverging  part 
of  the  nozzle  is  modeled,  and  its  geometry  approximates  one  of  the  two  profiles  used  in  Ref.  11.  The  nozzle 
throat  diameter  is  2.5  mm,  the  length  of  the  diverging  part  is  50.7  mm,  and  the  half-angle  is  20deg.  The 
surface  temperature  of  the  nozzle  is  assumed  to  be  300  K.  Numerical  results  are  obtained  for  two  throat- 
diameter  based  Reynolds  numbers,  1,230  (Case  I)  and  12,300  (Case  II),  with  the  stagnation  temperature 
of  333  K.  In  all  numerical  approaches,  the  computational  domain  starts  at  the  nozzle  throat  and  covers 
the  entire  diverging  part  of  the  nozzle,  as  well  as  a  small  part  of  the  plume  to  avoid  the  influence  of  the 
downstream  boundary  conditions.12  Four  different  numerical  approaches  are  used,  as  described  in  the  next 
section. 


III.  Computational  methods 


A.  DSMC  method 

The  DSMC-based  SMILE  computational  solver  was  used  in  this  work.  Details  on  the  tool  may  be  found 
elsewhere.13  The  SMILE  tools  that  were  used  in  the  present  work  include  axisymmetric  capability  with  radial 
weights,  different  grids  for  collisions  and  macroparameters,  both  of  which  are  two-level  Cartesian  grids  that 
automatically  adapt  to  flow  gradients,  and  parallel  implementation  with  efficient  load  balancing  techniques. 
The  majorant  frequency  scheme  was  employed  for  modeling  molecular  collisions.  The  VHS  model  was  used 
for  modeling  molecular  interactions.  Diffuse  reflection  with  full  energy  accommodation  was  assumed  on  the 
nozzle  wall.  For  both  selected  Reynolds  numbers,  solutions  independent  of  grid,  time  step,  and  number  of 
particles  were  obtained.  For  Re=l,230,  there  was  virtually  no  difference  observed  between  solutions  obtained 
for  0.3  million  cells  with  1  million  molecules,  and  3  million  cells  with  10  million  molecules.  For  Re=12,300, 
this  was  true  for  numerical  parameters  up  to  an  order  of  magnitude  larger. 

B.  Solution  of  NS  equations 

A  commercial  code,  CFD+- 1-  has  been  used  in  this  work  to  solve  the  Navier-Stokes  equations.  CFD+- 1-14  is  a 
flexible  computational  fluid  dynamics  software  suite  for  the  solution  of  steady  and  unsteady,  compressible  and 
incompressible  Navier-Stokes  equations,  including  multi-species  capability  for  perfect  and  reacting  gases.  In 
this  work,  a  perfect-gas  compressible  Navier-Stokes  solver  was  used  with  second  order  spatial  discretization 
and  implicit  time  integration.  Second  order  velocity  slip  and  temperature  jump  conditions  were  imposed 
on  the  nozzle  wall.  A  supersonic  inflow  with  prescribed  parameters  were  applied  at  the  nozzle  throat,  and 
backpressure  of  1  Pa  was  imposed  at  the  outflow  boundaries.  A  symmetry  condition  was  defined  at  the 
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nozzle  axis.  The  results  presented  below  (both  flow  fields  and  computational  requirements)  were  obtained 
for  a  multi-block  rectangular  grid  with  a  total  of  14,400  nodes.  The  computations  were  also  conducted  for 
four  times  smaller  and  four  times  larger  numbers  of  nodes,  and  found  fully  grid-resolved  with  14,400  nodes. 

C.  Finite  Volume  method  for  BGK  and  ES-BGK  equations 

A  finite  volume  solver  SMOKE  developed  at  ERC  has  been  used  to  deterministically  solve  the  BGK  and 
ES-BGK  equations.  SMOKE  is  a  parallel  code  based  on  conservative  numerical  schemes  developed  by  L. 
Mieussens.15  A  second  order  spatial  discretization  with  axial  symmetry  is  used  along  with  implicit  time 
integration.  A  supersonic  inflow  condition  is  used  at  the  nozzle  throat,  and  vacuum  outflow  conditions  are 
set  at  the  outer  boundaries.  Fully  diffuse  reflection  with  complete  energy  accommodation  is  applied  at  the 
nozzle  surface.  The  spatial  grid  convergence  was  achieved  increasing  the  number  of  nodes  from  3,600  to 
over  17,000.  The  convergence  on  the  velocity  grid  is  also  studied,  with  the  number  of  (x,  r,  0)  points  ranging 
from  (20,10,18)  to  (30,35,50). 

D.  Statistical  method  for  BGK  equation 

In  our  earlier  work,16  we  have  developed  and  studied  a  statistical  technique,  called  eDSMC,  which  models 
continuum  flows  through  a  collision  enforcement  procedure  which  guarantees  full  relaxation  of  the  molecular 
thermal  velocities  to  the  state  of  local  equilibrium.  The  technique  is  able  to  solve  inviscid  flows  with  tangency 
boundary  conditions  at  the  wall,  but  tends  to  under  predict  the  viscous  effects  in  the  boundary  layer  when 
diffusional  boundary  conditions  are  used.  In  this  work  we  continue  our  effort  to  apply  statistical  methods  to 
moderate  and  high  Reynolds  number  nozzle  flows  by  making  use  of  of  the  BGK  17  model.  It  is  interesting 
to  determine  whether  the  usual  DSMC  method  numerical  requirements,  that  the  cell  size  should  be  smaller 
than  the  local  mean  free  path  and  the  time  step  should  be  smaller  than  the  mean  collision  time,  could  be 
substituted  with  standard  CFD  criteria  such  as  the  Currant  number  and  the  resolution  of  flow  gradients. 

Recently,  a  number  of  authors8,10’18  have  developed  particle  approaches  to  the  solution  of  the  BGK 
and  ES-BGK  model  equations.  While  these  approaches  differ  in  details  and  theoretical  arguments  about 
the  formal  connection  between  the  numerical  schemes  and  the  equations  continue,  the  basic  ideas  of  the 
numerical  schemes  remain  the  same.  The  essence  of  these  kinetic  approches  is  to  select  a  given  number 
of  simulated  particles  from  those  available  in  a  computational  cell,  and  to  assign  to  them  new  velocities 
according  to  the  local  Maxwellian  (or  ellipsoidal)  distribution  function.  If  the  collision  frequency  for  such  a 
velocity  reassignment  is  properly  computed  and  local  values  of  the  translational  temperature  in  the  cell  are 
known,  then  the  procedure  certainly  mimics  the  collision  term  on  the  right  hand  side  of  the  BGK  equation, 

d(nf)/dt  +  c  •  d(nf)/d r  =  i m(fe  -  /) 

where  n  is  the  number  density,  v  is  the  collision  frequency,  /  is  the  molecular  velocity  distribution  function, 
fe  is  the  Maxwellian  distribution  function,  c  is  the  velocity  vector,  and  r  is  the  position  vector. 
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There  are  arguments  about  the  importance  of  certain  features  of  the  algorithms  to  preserve  energy  and 
momentum  conservation,  which  in  the  opinion  of  the  authors  is  a  secondary  problem,  since  on  average  these 
quantities  are  preserved.  The  major  issue  seems  to  be  the  absence  of  the  formal  derivation  of  the  numerical 
schemes  from  the  BGK  equation,  in  a  similar  manner  as  it  was  conducted  for  DSMC  schemes.19  In  this 
work  we  compare  results  of  a  statistical  BGK  scheme  with  those  of  the  finite  volume  solution  to  the  same 
equation.  We  believe  that  such  a  comparison  with  the  results  of  the  formally  derived  computational  scheme 
can  provide  some  insight  into  the  features  of  the  particle  BGK  scheme  which  we  use,  and  probably  to  particle 
BGK  schemes  in  general. 

The  details  of  the  statistical  BGK  scheme  are  as  follows.  The  collision  frequency  is  calculated  as 

v  =  Pr  -nk  T“_1 

V  Mre/  ) 

where  Pr  is  the  Prandtl  number  (1  for  the  BGK  equation,  in  this  work  we  used  both  1  and  0.7  with  no  visible 
difference  in  the  results),  k  is  the  Boltzmann  constant,  T  is  the  translational  temperature  in  the  flow,  /zre/  is 
the  gas  dynamic  viscosity  at  Tref.  The  collision  frequency  is  calculated  for  each  computational  cell  at  each 
time  step  based  on  the  local  translational  temperature  T  and  the  local  number  density  n.  The  local  number 
density  n  was  averaged  over  a  large  number  of  computational  time  steps,  while  the  local  temperature  T  is 
computed  based  on  the  instantaneous  velocities  of  the  computational  particles  in  the  cell. 

The  number  of  particles  in  a  cell  preselected  for  velocity  re-sampling  was  calculated  as  follows: 

Nc  =  int(AT(l  —  exp(— isAt)) 

where  At  is  the  time  step,  N  is  the  number  of  particles  in  the  cell  and  int  operator  means  the  nearest  smaller 
integer.  To  compensate  for  the  systematic  error  which  such  an  operator  produces,  one  more  particle  was 
added  to  to  the  list  of  preselected  particles  with  the  probability 

Pc  =  N(1  —  exp(— vAt)  —  int(iV(l  —  exp(— vAt)) 

The  preselected  particles  received  new  velocities  according  to  a  Maxwellian  distribution  at  the  local  cell 
temperature  and  velocity.  The  velocities  of  particles  which  have  not  been  preselected  remained  unchanged. 

Although  the  technique  is  not  limited  to  the  simple  gases,8  argon  was  used  as  the  working  gas  in  or¬ 
der  to  understand  the  basic  features  and  limitations  of  the  method  without  the  additional  complication  of 
translational-rotational  relaxation.  In  order  to  eliminate  the  impact  of  the  number  of  particles  per  computa¬ 
tional  cell,  a  large  number  of  particles  was  used  so  that  the  minimum  number  of  particles  per  cell  was  about 
150. 


IV.  Results  and  discussion 


A.  Low  pressure  Case  I 

The  results  presented  below  for  the  low  pressure  Case  I  were  obtained  by  a  DSMC  code,  its  modification 
implementing  the  statistical  BGK  scheme,  a  finite  volume  (FV)  ES-BGK  solver,  and  a  Navier-Stokes  solver. 
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Let  us  first  compare  the  results  obtained  by  FV  ES-BGK.  Again,  the  DSMC  result  is  assumed  to  the  the 
benchmark  result,  since  the  flow  condition  allowed  us  to  strictly  satisfy  all  DSMC  requirements,  and  the 
independence  of  the  DSMC  results  on  the  parameters  of  the  approach  was  verified.  Figures  f  and  2  show 
contour  plots  of  the  translational  temperature  and  axial  velocity  component,  respectively,  obtained  by  the 
two  numerical  schemes.  While  the  general  flow  structure  is  very  similar,  the  FV  BGK  scheme  was  found 
to  slightly  overpredict  the  viscous  effects  in  the  rarefied  portion  of  the  nozzle  flow.  This  is  seen  both  in 
temperature  and,  to  somewhat  larger  extent,  in  velocity.  Still,  the  total  difference  in  the  results  obtain  by 
different  approaches  is  very  small,  and  does  not  exceed  about  3%  in  any  portion  of  the  flow.  One  possible 
reason  for  the  difference  may  be  due  the  difference  in  the  collision  term  of  the  model  kinetic  equation  and 
the  full  Boltzmann  equation.  The  BGK  and  ES-BGK  equations  are  obtained  with  the  assumption  that  the 
collision  frequency  does  not  depend  on  velocities  of  colliding  particles,  and  that  the  particles  come  to  local 
equilibrium  (or  to  the  ellipsoidal  distribution)  after  they  experience  vtc  collisions,  where  v  and  rc  are  the 
collision  frequency  and  collision  time,  respectively.  Note  that  in  addition  to  ES-BGK,  the  BGK  equation 
was  also  solved  using  a  finite  volume  method  for  Case  I,  and  the  difference  in  the  results  between  the  BGK 
and  ES-BGK  solutions  found  to  be  within  a  percent. 

After  having  established  the  accuracy  of  the  FV  ES-BGK  scheme,  let  us  now  compare  the  results  obtained 
with  different  approaches  to  the  solution  of  model  kinetic  equations.  Figures  3  and  4  show  the  contour 
plots  of  the  translational  temperature  and  axial  component  of  velocity  obtained  by  the  deterministic  and 
statistical  approaches.  The  agreement  is  rather  good  in  the  first  quarter  of  the  nozzle,  while  in  the  remaining 
part  the  particle  BGK  tends  to  overpredict  the  viscous  effects.  For  the  particle  solution,  the  axial  velocity 
near  the  nozzle  plane  is  lower,  and  the  temperature  is  higher,  than  for  the  finite  volume  solution.  The 
difference  between  the  particle  BGK  and  DSMC  is  even  greater.  Note  that  similar  effect  was  observed  in 
Ref.  8  for  a  diatomic  gas.  Possible  reasons  for  such  a  difference  require  further  investigation. 

The  comparison  of  the  FV  ES-BGK  results  with  the  NS  predictions  is  given  in  Figs.  5  and  6.  While  the 
temperature  contours  are  in  good  agreement  near  the  nozzle  centerline,  there  is  a  large  difference  both  in 
translational  temperature  and  axial  velocity  near  the  wall.  The  NS  solution  is  characterized  by  significantly 
lower  velocity  and  higher  temperatures  in  that  region,  with  the  difference  increasing  toward  the  nozzle  exit 
plane.  The  main  reason  for  such  a  difference  is  that  the  velocity  slip  and  temperature  jump  boundary 
condition  used  in  the  NS  solver  does  not  provide  a  correct  description  of  gas-surface  interaction,  and  the 
difference  increases  with  the  degree  of  thermal  nonequilibrium. 

An  illustration  of  thermal  nonequilibrium  in  the  flow  is  provided  in  Fig.  7,  where  the  ratio  of  the  tem¬ 
perature  in  axial  direction  to  the  overall  temperature  is  presented.  In  equilibrium,  this  ratio  should  be  equal 
to  the  unity,  and  the  difference  from  the  unity  serves  as  the  measure  of  the  flow  thermal  nonequilibrium. 
The  temperature  ratio  presented  in  this  figure  was  obtained  with  the  DSMC  method.  The  results  show  that 
there  is  a  strong  nonequilibrium  near  the  nozzle  lip,  where  the  ratio  is  as  low  as  0.7.  Generally,  a  noticeable 
degree  of  nonequilibrium  exists  near  the  nozzle  wall,  where  the  axial  temperature  is  lower  than  the  radial 
temperature,  and  in  the  middle  of  the  nozzle  off  the  centerline,  where  the  radial  temperature  becomes  larger. 
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The  flow  nonequilibrium  and  the  difference  between  the  translational  mode  temperatures  indicate  that  the 
accuracy  of  the  NS  equations  becomes  questionable,  and  the  may  no  longer  be  applicable.  Comparison  of 
the  NS  and  DSMC  temperature  contours  given  in  Fig.  8  clearly  shows  the  impact  of  the  flow  nonequilibrium: 
there  is  virtually  no  difference  between  the  two  solutions  near  the  nozzle  axis,  whereas  near  the  nozzle  wall 
they  significantly  differ. 

To  quantify  the  differences  between  different  approaches,  the  gas  parameters  are  presented  along  two 
cross  sections,  the  nozzle  centerline  and  the  nozzle  exit  plane.  The  density,  velocity,  and  temperature 
along  the  centerline  are  shown  in  Figs.  9-11.  Here,  X=0  corresponds  to  the  nozzle  throat.  One  of  the 
conclusions  here  is  that  three  of  four  our  approaches  produce  very  similar  results  for  all  gas  properties  under 
consideration.  Another  conclusion  is  that  the  fourth  approach,  namely  particle  BGK,  produces  almost  4% 
lower  flow  velocites  and  25%  larger  gas  temperatures  closer  to  the  nozzle  exit  plane.  The  reasons  for  the 
deviation  needs  to  be  investigated,  here  we  note  only  that  the  impact  of  the  number  of  particles,  time  step, 
or  the  cell  size  in  the  particle  BGK  method  is  unlikely  to  be  noticeable,  since  the  convergence  of  results  to 
these  parameters  has  been  tested. 

The  axial  velocity  and  temperature  along  the  nozzle  exit  plane  are  presented  in  Figs.  12  and  13.  Here, 
Y=0  corresponds  to  the  nozzle  centerline.  The  VF  ES-BGK  solution  is  close  to  that  of  DSMC  both  for  the 
axial  velocity  and  gas  temperature.  The  particle  BGK  result  agrees  with  the  DSMC  one  near  the  wall,  but 
start  to  deviate  noticeably  for  smaller  radial  coordinates.  The  NS  solver  produces  results  accurate  in  the 
coreflow  but  strongly  different  from  the  ES-BGK,  DSMC  and  particle  BGK  throughout  the  boundary  layer. 
Generally,  we  can  conclude  that  the  FV  ES-BGK  results  are  in  very  good  agreement  with  the  DSMC  in  the 
entire  diverging  part  of  the  nozzle.  The  solution  of  the  model  kinetic  equation  accurately  captures  both  the 
boundary  layer  and  the  nozzle  coreflow.  The  particle  BGK  results  deviate  from  DSMC  in  the  coreflow,  and 
the  Navier-Stokes  equations  appear  to  be  inaccurate  in  the  boundary  layer. 

The  second  most  important  consideration  for  the  selection  of  the  appropriate  numerical  technique  for 
modeling  nozzle  flows,  in  addition  to  physical  accuracy,  is  computational  efficiency.  The  summary  of  the 
computational  requirements  necessary  to  complete  the  simulations  described  in  this  section  are  shown  in 
Table  1.  Note  that  while  the  CPU  time  given  in  this  table  represents  the  actual  computational  cost  of  the 
runs,  some  caution  need  to  be  used  when  analyzing  the  results.  The  particle  simulations  where  conducted 
until  the  statistical  scatter  in  major  macroparameters  was  not  visible  inside  the  nozzle.  This  typically 
required  about  10,000  timesteps  to  reach  steady  state  and  over  100,000  more  time  steps  for  sampling.  It  is  of 
course  possible  to  significantly  reduce  the  number  of  sampling  steps  if  larger  statistical  scatter  is  acceptable. 
Decreasing  the  number  of  particles  is  possible,  but  will  not  decrease  the  computational  time,  since  a  larger 
number  of  time  steps  for  sampling  will  be  needed.  The  finite  volume  BGK  run  time  may  be  decreased  by  up 
to  a  factor  of  two  without  significant  loss  in  accuracy,  since  it  is  possible  to  use  somewhat  smaller  numbers 
of  spatial  cells  and  velocity  grid  points. 

The  results  show  that  the  Navier-Stokes  solver  is  at  least  two  orders  of  magnitude  more  efficient  than  the 
kinetic  solvers,  although  its  accuracy  is  questionable,  as  was  illustrated  earlier  in  this  section.  The  DSMC 
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method  is  faster  than  the  other  two  kinetic  approaches,  so  for  the  problem  under  consideration  the  use 
of  a  BGK  code  does  not  seem  to  be  justified.  Note  also  that  the  DSMC  method  may  requires  somewhat 
less  intervention  from  the  user,  since  satisfying  the  usual  DSMC  requirements  on  the  time  step,  number  of 
particles,  and  cell  size  provides  accurate  results  without  additional  convergence  study.  For  FV  ES-BGK, 
a  convergence  analysis  for  cell  size  and  velocity  space  may  be  necessary  unless  excessively  fine  spatial  and 
velocity  grids  are  used. 

B.  High  pressure  Case  II 

Analysis  of  Case  I  has  showed  that  the  FV  ES-BGK  is  comparable  to  the  DSMC  method,  both  in  accuracy. 
The  solution  of  the  NS  equations,  while  orders  of  magnitude  more  efficient,  strongly  underestimates  flow 
velocity  and  overestimates  temperature  in  the  boundary  layer.  The  NS  equations  are  expected  to  be  much 
more  accurate  for  Case  II,  where  the  density  is  an  order  of  magnitude  larger  than  in  Case  I.  At  the  same  time 
the  application  of  the  DSMC  method  to  Case  II  is  much  more  time  consuming,  and  the  conventional  DSMC 
requirements  were  not  strictly  satisfied  even  when  over  50  million  particles  and  15  million  cells  were  used. 
The  comparison  of  the  DSMC  results  obtained  for  smaller  number  of  particles  and  cells  and  coarser  time 
steps  have  shown  however  that  the  above  DSMC  parameters  are  sufficient  for  obtaining  numerical  parameter 
independent  DSMC  results.  Significant  computational  difficulties  were  encountered  when  applying  the  FV 
ES-BGK  to  this  higher  density  case,  as  several  times  larger  number  of  velocity  points  are  required  to  obtain 
grid  converged  results.  Details  about  the  accuracy  and  efficiency  of  the  above  numerical  approaches  are 
presented  below. 

In  addition  to  the  four  approaches  used  to  calculate  Case  I,  the  eDSMC  technique16  was  also  used  for  Case 
II.  The  collision  enforcement  eDSMC  technique  makes  the  flow  effectively  inviscid.  It  potentially  requires 
less  computational  cells  and  particles  and  larger  time  steps  than  both  DSMC  and  statistical  BGK  methods. 
The  time  step  in  the  eDSMC  technique  does  not  depend  on  the  local  mean  free  path  but  is  rather  defined  by 
the  CFD  type  CFL  number.  When  a  boundary  condition  of  specular  reflection  of  computational  molecules 
from  the  walls  is  used,  eDSMC  is  able  to  properly  solve  the  inviscid  nozzle  flows.  In  Ref.16  it  was  shown 
that  the  efficiency  of  the  baseline  DSMC  method  could  be  improved  when  a  combination  of  eDSMC  and 
DSMC  is  applied.  Although  the  technique  underestimates  the  viscous  effects  in  the  boundary  layer,  when 
the  diffusional  boundary  conditions  are  used,  it  is  still  capable  of  adequately  solving  the  core  portion  of  a 
nozzle  flow. 

A  side-by-side  comparison  of  the  flowfields  for  Case  II  is  given  in  Figs.  14-19.  Similar  to  Case  I,  the  FV  ES- 
BGK  solution  is  in  fairly  good  agreement  with  the  DSMC  prediction  throughout  the  computational  domain. 
The  particle  BGK  results  predict  velocity  and  temperature  values  similar  to  BGK  in  the  downstream  portion 
of  the  coreflow.  However,  the  boundary  layer  is  visibly  thicker,  and  a  weak  compression  front  coming  from 
the  wall  at  the  nozzle  throat  to  the  centerline  at  approximately  6  mm  from  the  throat,  predicted  by  both 
DSMC  and  FV  ES-BGK,  is  not  captured  by  particle  BGK.  The  most  probably  reason  for  these  problems  is 
insufficient  cell  resolution  near  the  nozzle  wall.  The  solution  of  the  NS  equation  is  rather  close  to  the  FV 
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ES-BGK-one,  with  a  relatively  small  difference  observed  in  the  boundary  layer. 

A  quantitative  comparison  of  flow  velocity  and  temperature  for  Case  II  is  given  in  Figs.  20  and  21,  where 
the  profiles  along  the  centerline  are  shown,  and  in  Figs.  22  and  23  for  the  cross  section  along  the  nozzle 
exit  plane.  It  is  clearly  seen  that  the  DSMC,  NS,  and  ES-BGK  profiles  of  both  velocity  and  temperature 
are  very  close  along  the  nozzle  centerline.  The  particle  BGK  solution  deviates  from  the  other  three  in  the 
region  where  the  compression  wave  comes  to  the  centerline,  and  and  agrees  with  them  otherwise.  The  results 
obtained  with  eDSMC  are  similar  to  particle  BGK  in  the  coreflow.  The  results  along  the  exit  plane  show 
that  there  is  still  some  difference  between  the  NS  and  DSMC  predictions,  both  in  temperature  and  flow 
velocity.  As  in  Case  I,  this  difference  is  explained  by  the  limitations  of  the  velocity  slip  boundary  condition 
used  in  NS.  There  is  also  some  impact  of  flow  nonequilibrium,  but  its  contribution  is  significant  only  near 
the  nozzle  lip,  where  Tx/T  ration  reaches  0.85.  Figures  22  and  23  show  that  the  eDSMC  results  are  close 
to  the  results  of  DSMC  and  other  methods  in  the  core  of  the  flow,  where  viscosity  effect  are  not  significant, 
but  deviate  greatly  from  DSMC  in  the  boundary  layer  region. 

The  results  that  summarize  the  numerical  efficiency  of  different  approaches  are  shown  in  Table  2.  The 
deterministic  solution  of  the  model  kinetic  equations  becomes  noticeably  less  efficient  than  the  other  two 
kinetic  approaches,  since  a  much  larger  number  of  spatial  cells  and  velocity  grid  point  is  needed  to  obtain  a 
converged  result.  As  much  as  50  angular  velocity  cells  are  necessary  in  this  case,  compared  to  less  than  20  for 
Case  I.  It  is  interesting  to  note  that  computational  efforts  for  the  DSMC  and  statistical  BGK  schemes  were 
comparable,  but  the  convergence  process  was  quite  different.  The  DSMC  schemes  tend  to  converge  faster, 
but  requires  more  computational  effort  to  collect  the  sufficient  information  for  the  solution  to  be  smooth. 
Slower  convergence  of  the  statistical  BGK  method  is  probably  due  to  the  ’’history”  of  the  macroparameter 
sampling  procedure  which  defines  the  collision  frequencies.  This  history  may  at  the  same  time  explain  smaller 
statistical  scatter  of  particle  BGK.  Similar  to  Case  I,  the  NS  solver  was  over  two  order  of  magnitude  more 
efficient  than  the  kinetic  methods,  however,  the  NS  results  still  differ  from  DSMC.  The  eDSMC  method 
was  an  order  of  magnitude  faster  than  the  other  statistical  methods,  which  makes  it  a  promising  statistical 
approach  for  solving  the  core  portion  of  the  nozzle  flow  in  a  combination  with  other  statistical  methods  to 
proper  resolve  the  boundary  layer. 


V.  Conclusions 

Argon  flow  through  a  conical  nozzle  was  studied  for  two  Reynolds  number,  1,230  and  12,300,  using  four 
different  approaches.  These  include  one  continuum  approach  (solution  of  Navier-Stokes  equations)  and  three 
kinetic  approaches  (the  DSMC  method  and  a  statistical  and  deterministic  methods  for  the  model  kinetic 
equations).  Analyses  of  the  accuracy  of  the  approaches  and  their  numerical  efficiency  was  conducted.  Several 
conclusion  can  be  drawn  from  the  results  of  the  computations.  The  most  accurate  results  in  both  high  and 
low  pressure  cases  are  believed  to  be  those  obtained  the  DSMC  method.  All  standard  DSMC  requirements 
were  satisfied  for  Re=l,230.  For  the  higher  pressure  case,  although  parameters  of  the  numerical  scheme  were 
somewhat  relaxed,  a  solution  independent  on  the  time  step,  number  of  simulated  molecules,  and  cell  size, 
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was  obtained. 

The  Navier-Stokes  solutions  are  in  good  agreement  with  the  DSMC  results  in  the  higher  density  portion 
of  the  flow  and  in  the  coreflow,  where  rarefaction  effects  are  small.  In  the  boundary  layer,  even  though  the 
velocity  slip  and  temperature  jump  boundary  conditions  were  used,  there  is  a  difference  between  the  NS  and 
DSMC  solutions.  The  finite  volume  solution  of  the  ES-BGK  equation  is  in  very  good  agreement  with  the 
DSMC  method  in  the  entire  computational  domain  for  both  Reynolds  numbers.  The  particle  BGK  method 
agrees  with  the  DSMC  results  in  the  boundary  layer  for  the  lower  pressure  case  and  in  the  coreflow  for  the 
higher  pressure  case.  Although  further  analysis  is  needed,  the  deviation  of  the  statistical  BGK  results  from 
DSMC  is  currently  attributed  to  the  lack  of  formal  definition  of  the  particle  BGK  scheme,  which  may  result 
in  solutions  different  from  the  exact  solution  of  the  BGK  model  equation.  Some  of  the  problems  in  the  high 
pressure  areas  of  the  flow  can  be  attributed  to  the  grid  resolution,  although,  it  was  observed  to  be  adequate 
in  the  rarefied  portion  of  the  flow  where  the  largest  difference  is  seen. 
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Table  1.  Parameters  of  the  methods,  Case  I 


Method 

DSMC 

Part  BGK 

FV  ES-BGK 

NS 

Time  (CPUH) 

200 

1000 

400 

<1 

Number  of  particles 

10  mil 

50  mil 

- 

- 

Number  of  cells 

3  mil 

0.5  mil 

3,600 

3,600 

Table  2.  Parameters  of  the  methods,  Case  II 


Method 

DSMC 

Part  BGK 

FV  ES-BGK 

NS 

eDSMC 

Time  (CPUH) 

1000 

1000 

5000 

3 

100 

Number  of  particles 

50  mil 

50  mil 

- 

- 

1  mil 

Number  of  cells 

15  mil 

0.5  mil 

17,200 

14,400 

100,000 
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Figure  1.  Case  I,  DSMC  and  FV  ES-BGK  temperature  contours  (K). 


Figure  2.  Case  I,  DSMC  and  FV  ES-BGK  axial  velocity  contours  (m/s). 
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Figure  3.  Case  I,  Particle  BGK  and  FV  ES-BGK  temperature  contours  (K). 


Figure  4.  Case  I,  Particle  BGK  and  FV  ES-BGK  axial  velocity  contours  (m/s). 
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Figure  5.  Case  I,  NS  and  FV  ES-BGK  temperature  contours  (K). 


Figure  6.  Case  I,  NS  and  FV  ES-BGK  axial  velocity  contours  (m/s). 
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0.015 


Figure  7.  Case  I,  DSMC  Tx/T  ratio. 


Figure  8.  Case  I,  DSMC  and  NS  temperature  contours  (K). 
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X  [m] 


Figure  9.  Case  I,  Density  along  the  centerline. 


Figure  10.  Case  I,  Axial  velocity  along  the  centerline. 


17  of  24 


American  Institute  of  Aeronautics  and  Astronautics 
Distribution  A:  Approved  for  public  release;  distribution  unlimited 


X  [m] 


Figure  11.  Case  I,  Temperature  along  the  centerline. 
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Figure  12.  Case  I,  Axial  velocity  at  the  nozzle  exit. 


T 


Figure  13.  Case  I,  Temperature  at  the  nozzle  exit. 
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Figure  14.  Case  II,  DSMC  and  FV  ES-BGK  temperature  contours  (K). 
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Figure  15.  Case  II,  DSMC  and  FV  BGK  axial  velocity  contours  (m/s). 
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Figure  16.  Case  II,  Particle  BGK  FV  ES-BGK  temperature  contours  (K). 


X 

Figure  17.  Case  II,  Paricle  and  FV  BGK  axial  velocity  (m/s). 
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Figure  18.  Case  II,  NS  and  FV  ES-BGK  temperature  contours  (K). 


X 


Figure  19.  Case  II,  NS  and  FV  ES-BGK  axial  velocity  contours  (m/s). 
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Figure  20.  Case  II,  Axial  velocity  along  the  centerline. 


X  [m] 


Figure  21.  Case  II,  Temperature  along  the  centerline. 
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Figure  22.  Case  II,  Axial  velocity  at  the  nozzle  exit. 


Figure  23.  Case  II,  Temperature  at  the  nozzle  exit. 
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